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Abstract 

An accurate and systematic equation of state for the hard-core one-component 
plasma (HCOCP) is obtained. The result is based on the Hubbard- Schofield trans- 
formation which yields the field-theoretical Hamiltonian, with coefficients expressed 
in terms of equilibrium correlation functions of the reference hard-core fluid. Explicit 
calculations were performed using the Gaussian approximation for the effective Hamil- 
tonian and known thermodynamic and structural properties of the reference hard-core 
fluid. For small values of the plasma parameter T and packing fraction the Debye- 
Hiickel result is recovered, while for F S> 1, the excess free energy Fex Sind internal U^x 
energy depend linearly on F. The obtained expression for Uex is in a good agreement 
with the available Monte Carlo data for the HCOCP. We also analyse the validity of the 
widely used approximation, which represents the free energy as a sum of the hard-core 
and electrostatic part. 



Preprint accepted to Journal of Physics D: Applied Physics. 



1 Introduction 

The one component plasma (OCP) is one of the basic models in the field of charged systems. 
The OCP model is formulated as a system of point particles, interacting via the Coulomb 
potential, which move in a uniform neutralizing background. It has found important phys- 
ical applications in a variety of fields ranging from terrestrial physics, through important 
technological applications to cosmology [|I], ^ ^ || . As a reference model it is used in many 
areas of soft condensed matter, such as colloidal and polyelectrolyte solutions, e.g. p, |^, |^, |^, 
etc. All thermodynamic properties of the OCP depend only on the dimensionless plasma 
parameter F = Ib/o-c, where Ib = e^/kBT is the Bjerrum length (e is the charge of the 
particles, /cb is the Boltzmann's constant, T is the temperature) and Oc = (3/47rp)^/^ is the 
ion-sphere radius with p = N/Q being the concentration of particles (A^ is the number of 
particles, fl is the volume of the system). Using a field-theoretical approach, a fairly accurate 
and simple expression for the equation of state of the OCP has been obtained within the 
Gaussian approximation for the effective Hamiltonian contrary to previous calculations, 
e.g. pO] , pT| ], this gives a correct behavior for the thermodynamic functions in the full range 



of r (see and references therein) and does not have fitting parameters as e.g. in [11, 12]. 
In |]13[ field-theoretical calculations for the equation of state, going beyond the Gaussian 
approximation, have been performed, showing that corrections to the Gaussian theory are 
rather small. 

A closely related model - the hard-core one component plasma (HCOCP) incorporating 
a hard core repulsion between ions - gives a more satisfactory description of the short-range 
electrostatic correlations. The importance of this model follows also from the fact that it 
belongs to the class of the so-called primitive models used to describe molten salts [Q, 
electrolytes |]15[, liquid metals [|1^, |1^, |18| and charged colloidal solutions [|1^, The 
mean spherical approximation (MSA) was applied to account for both the hard-core and 
electrostatic interactions in this system p3|, |2 



Although it has an analytical solution, it 
exhibits sometimes unphysical negative contact values for the pair correlation function, which 
has to be remedied by a rescaling procedure A mixed Percus-Yevick/hypernetted 

chain integral equation also has been used pSl and it was observed that it fits better the 



simulation results than the MSA and cluster expansion |31|]. The former however 
does not have an analytical formulation. An analytical equation of state for the HCOCP 
has been proposed in [^, as a simple generalization of the hole-corrected Debye-Hiickel 
theory |10|, which in addition to the correlation hole around charged particles takes into 
account the hard-core repulsion. This was used afterwards to develop a generalized van 
der Waals theory p3i. The basic physical idea exploited in this theory is that due to the 



strong electrostatic repulsion a "hole" appears around a charged particle from which all other 
particles are expelled. Outside the hole the electrostatic interactions are not very strong and 
may be described on the Debye-Hiickel level; the size of the hole is found self-consistently. 
Although being physically appealing, this theory does not give a satisfactory description 
for large values of the plasma parameter P for the OCP |jlO[, and is also not accurate for 
the HCOCP for large packing fractions rj (for rj = 0.4 the deviations from MC data for the 
internal energy reach 24% 133]). Recently an exact low-density expansion for the free energy 



of the HCOCP has been obtained p4 
electrostatic interaction 



This however can not be applied in the case of strong 
i.e. for large P. 

In the present study we develop a theory which allows to derive a fairly accurate and 
simple equation of state for the HCOCP in the whole range of P from the Debye-Hiickel 
limit P -C 1 up to the limit of strong coupling P ^ 1. It reproduces within 1-3% accuracy 
the available Monte Carlo (MC) data for 0.2 < P < 70. Larger deviations occur in the 
region where the MC data are not very accurate (see the discussion below). As in the case 
of the OCP we use the Hubbard-Schofield transformation which yields a field-theoretical 
Hamiltonian for the HCOCP with coefficients expressed in terms of equilibrium correlation 
functions of a reference hard-core fluid. Although we perform the explicit calculations within 
the Gaussian approximation for the effective Hamiltonian, one can go beyond the Gaussian 
approximation using the standard perturbation technique, provided the structural properties 
of the reference hard-core fluid are known. Since the obtained equation of state is in good 
agreement with available MC data we use it to address the problem of additivity of the hard- 
core and electrostatic contributions to the excess free energy. This seems to be important 
since it is usually assumed that these contributions are additive (e.g. ||35|) and the range of 
validity of such approximation has not been studied yet. 

The rest of the paper is organized as follows: in Sec. 2 we briefly sketch the Hubbard- 
Schofield approach which yields the field-theoretical Hamiltonian for the HCOCP. In Sec. 3 we 
derive the equation of state within the Gaussian approximation for the effective Hamiltonian, 
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and compare the theoretical findings with available Monte Carlo data. In that section we 
also analyse the accuracy of the widely used approximation in which the excess free energy 
is written as a sum of the hard-core and electrostatic part. In the Conclusion we summarize 
our findings. 



2 Field-Theoretical Hamiltonian for HCOCP 

We start from the HCOCP Hamiltonian with omitted ideal part which may be written as 
follows ip-^ = knT): 

H = \p-^ J2' ^kip^P-^ -p) + Hue (1) 

where the first term in the right-hand side of (|I]) refers to the Coulomb interactions, written 
in terms of collective density variables, 

where denotes coordinate of j-th particle, 

z/fc = A-KlB/k^ (3) 

is the Fourier-transformed Coulomb potential and Hhc describes the hard-sphere interaction. 
Summation in (|I]) is to be performed over the wave vectors k = {kxkykz} with ki = 2'7ili/L 
{i = x,y,z), where k are integers, = fl, and the prime over the sum denotes that the 
term with k = is excluded 



2.1 Hubbard- Schofield transformation 



The configurational integral may be written in terms of the configurational integral of the 
reference (hard-sphere hquid) system Qji [R^, RS, 



as: 



g = ^exp |-i^'z/fc(pkp-k - Q 



(4) 



where ((. . .))^ = Q^^^ J dr^{. . .) denotes the averaging over the reference system. In accor- 
dance with the Hubbard-Schofield scheme we use the identity: 



exp [ 2'^^'''^ 



V2 



1 1 

/ exp ( — y"^ / + ixy ] dy 



and arrive after some algebra at: 



Q = QrJ n'ckC?<^kexp|-^^'z/fc VkV^-kj ^exp ji^'pkV'- 



(5) 



R 



where Ck = (27rz/fc)~^/^ exp{z/fcp/2}, and where the integration is to be performed under 
restriction, yj.k = v^k (v'k is the complex conjugate of Lpi^) [^. Applying the cumulant 
theorem [^] to the factor (exp {z X^'k PkV-k})/? one obtains: 
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Q = Qr I n''^k'iv5ke ^, with 

n=2 ki,...k„ 
M2(kl,k2) = ^5ki+k2,o|— + (PkiP-ki)c 



n 



(6) 



Mn(kl, 



—I 



here (. . .)^^ denotes cumulant average |^ for the reference hard sphere fluid system. Note 
that (|^) gives the field-theoretical expression for the partition function with, being the 
Fourier components of the scalar field ^{r) The coefficients of the effective Hamiltonian 
(^) are expressed in terms of the correlation functions of the reference hard-core fluid. Using 
definitions of correlation functions of fluids H3] and definitions of the cumulant averages 



4l| , one can directly evaluate (pki • • • Pk„)c_R (and thus the coefficients M„(ki, . . . k„)). It 



is straightforward to show that (pki ■■■P]i„)cr be expressed in terms of the Fourier 
transforms of the correlation functions h2, h^, . . . , hn of the reference system, defined as 

/i2(ri,r2) = 5f2(ri,r2) - 1 (7) 

h{ri, r2, rg) = g-i{ri, r2, rg) - g2{ri, r2) - 5f2(ri, rg) - g2{r2, rg) + 2 , etc., (8) 

where gi{ri,..,ri) - are Z-particle correlation functions In particular the second-order 
cumulant {pkP^k)^j^ may be written as: 



(PkP-k) 



cR 



(9) 



and the Fourier transform of the function /i2(k) is related to the direct correlation function 
52(k) 

/i2(k) = g2(k)/[l-p52(k)], (10) 



for which an explicit analytical expression is known |44 



2.2 Gaussian approximation for the effective Hamiltonian 

Now we concentrate on the Gaussian part of the effective Hamiltonian, i.e. we skip all the 
terms with a power of the field larger than two. The accuracy of this approximation has been 



critically examined in Ref. |13| for the OOP without hard-core interactions by perturbatively 
calculating higher-order terms. It was found that higher order terms contribute very little 
to the free energy and that the Gaussian approximation is in fact excellent over the whole 
range of coupling parameters. One can therefore assume, and this is indeed borne out by 
our comparison with Monte-Carlo data below, that the Gaussian approximation should also 
be quite good for the present case. Using (0) and (P) we write for this case: 

Performing (Gaussian) integration over y^k we arrive after some algebra at 



\ -1 
P \ 



pc2(k) 



(12) 
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Taking the logarithm of the configuration integral one obtains the free energy, and since the 



analytical expression for C2(k) is available |^^, no additional approximation is, in principle 



required. This leads, however, to an expression which is to be evaluated numerically. We note 
that owing to the long-range nature of the Coulombic interactions, the main contribution to 
the free energy comes from the long-wave modes of the density fluctuations, which correspond 
to small k. Therefore only the small- A; behavior of the direct correlation function is important. 
This suggests to approximate C2(k) by a truncated expansion 

C2(k)~g2(0)-£2(0)"A;2 + ... , (13) 

which correctly behaves at small k. As we show in what follows, only wave-vectors with 
k < ko contribute to the configuration integral, thus we apply the approximation (|T^ for 
the interval < k < ko. By numerical evaluation of the free energy using the full expression 
for 62 (k), we convinced ourselves that deviations of the quadratic form from the actual C2(k) 
at larger k do not noticeably affect the results [^]. We therefore do not require k and ko to 
be small since the particular behavior of C2(k) for large k is not important. On the other 
hand, the quadratic approximation allows to obtain an analytical equation of state for the 
HCOCP that reproduces fairly well the available Monte Carlo data. With (0) one can write 
for the configurational integral: 

fc^>0 Vk 1-PC2(0) (l-pC2(0))2y 

The most accurate estimate for 62(0) may be found, using 62(0) = /i2(0)/(l + ph2{0)) from 
(p!0|) and the relation for the isothermal compressibility Xr^ = p{9Pr/ dp)j3 |Q, where Pr is 



the pressure of the reference system |l47|| : 

l + ph2{0) = pkBTxR = Zo (15) 

The value of Zq follows from the fairly accurate Carnahan-Starling free energy of hard-sphere 
fluid dSf: 

^=(1^^' ^''^ 

where 77 = ^ixpd? is the packing fraction and d is the diameter of the spheres; this yields Zq 
as a second derivative of F^c with respect to density: 

Zo = (1 - r])\l + 4:7] + V - 4r?^ + ^')"' (17) 

Using the Wertheim-Thiele solution for the direct correlation function C2(r) |Q and defini- 
tion of 02(0)" from Eq.(0) one obtains: 

^2(0)" = ^ y r'^C2{r)dT = -(7rc/Vl20)(16 - Ur] + V)(l - r/)~^ (18) 

With g"(0) from (|I|) and £2(0) = (Zo-l)/(pZo) from (^) and (|15|) we arrive at the following 
expression for the configurational integral 

Q = Qr n' exp(pz/fc) {pukZo + 9)-' (19) 



where 9 = 1- ATrlBipZofd'^iO). 
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3 Equation of state for the HCOCP 



3.1 Equation of state and its accuracy 



Now we show that within the Gaussian approximation to the effective Hamihonian, one 
obtains rather accurate equation of state for HCOCP provided that an appropriate value of 
the "ultraviolet cutoff" in the k-space is employed. From (|^) we find for the excess free 
energy of the HCOCP: 



- = ln(Q) 



-pF'hc + \ E'b^fc - Hp^kZQ + 6)] 



(20) 



Now we argue that the summation in ( ^DD should be carried out over a finite number of the 



wave-vectors k. In this we follow the Debye theory of the specific heat of solids (e.g. p9 



Namely, we assume that the total number of degrees of freedom in the system, 3iV, should 
be equal to the total number of physically different modes with the wave-vectors k within 
the spherical shell of radius fco in the k-space. The number of modes is twice the number 
of the wave-vectors, since for each k one has a sine and cosine mode (the amplitude of the 
k-th mode is a complex number) ||5^. Thus we obtain: 



2 — -A-n 



ko 



Pdk = 3N 



(21) 



where the factor Q/Sn^ appears when the integration in k-space is used instead of summation. 
From (l2l|) follows that 



A similar Debye- like scheme to find the cutoff fco was first proposed for plasma in ||5T[ , where 
a somewhat different value of the cutoff wave-vector was reported. Using as obtained 
above we write: 



N 



N 
N 



+ - 



2 87r3 N Jo 



+ 



4 Jo 



/ k dk [In (pUkZo + 0) - pz/fc] 
Jo 

bTZo\ 



In 



(22) 



where x = /c/Zcq, so that pz/^ = 

where we define the constant h 
for the excess free energy 



AtiIbp/ {k^xY = hV/x"^ (see (^) for the definition of z/^), and 

9 / 9 \ 

= g ( . The last integral is easily evaluated and yields 



ksTN 
with 



4?] - ?,if 



In (0 + fePZo) - 6r 3 



2Zo 




3 ( hTZo 



arctan 



= 1 + 



6 r/2(l -?7)4(16- llr/ + V 



5 dksT (1 + 4r/ + 4r/2 - Ar]^ + r]^y 
and for the excess internal energy, Uex = —T'^d{Fex/T)/dT of the HCOCP 




Uex 

ksTN 



bVZo 




br 



\ 



arctan 



/ 
bTZo 



(25) 
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For T] = 0, Eqs.(23.2^) recover the corresponding result for the one component plasma [Q. 

As it follows from ( p3D and (|25|), for F — > and 77 — > the Debye-Huckel behavior is 
obtained. In the opposite limit F ^ 1 and for any packing fraction rj Eqs.(^) and (|25| ) 
demonstrate a linear behavior on F. The leading term for this case is —AT, where the 
constant A reads: 

A = h = l(^y^' = 0.881... 



4 2 Vtt 



r2 



This is fairly close to the constant A = 0.899 ... of the fits for the OCP (see e.g. ||5^, |53|1). 

For arbitrary values of 77 and F we compare our analytical expression for the excess 
internal energy (|25| ) with the available MC data for the HCOCP (Fig.l). Fig. 2 shows the 
relative error of the analytical expression (|25|) . As it follows from Fig.l and Fig.2 the equation 
of state is fairly accurate in full range of plasma parameters for which MC data are available. 
For F > 1 the relative error does not exceed 1 — 3% for all values of packing fraction and 
plasma parameter. For F > 10 one observes the linear behavior and deviation from the 
numerical data for this range of F less than 1%. 

The maximal deviation of the analytical expression from the MC data occurs for 0.2 < 
F < 1. In this range of the plasma parameter the deviation is about 5% with the maximal 
one of 12% at the smallest value F = 0.2. It should be noted however, that such deviation 
occurs for F, where the method of MC loses its accuracy. Moreover, we expect that for very 
small packing fraction, t] = 0.001, where the maximal deviation is observed, the contribution 
to the internal energy due to the hard-core interactions may not exceed 1%. Therefore for 
this 7] the difference between Uex of the HCOCP and Uex of the OCP is less then 1%. It 
has been already shown that the analytical equation of state for the OCP (which follows 
from (^) for t] = 0) has the accuracy of 8% Thus we expect that the accuracy of our 
equation of state may not be worse than 9% for these values of 77; the observed deviation of 
12% seems to be the manifestation of the low accuracy of MC for this range of parameters. 



3.2 Additivity of the hard-core and electrostatic components of 
the free energy. 

With the equation of state obtained we can analyse the accuracy of the widely used ap- 
proximation for the excess free energy, which represents this as a sum of the hard-core and 



electrostatic component, F^c + Fex,ocp (see e.g. |]3^). The validity of this approximation 
may be checked by the direct comparison of the excess free energy of the HCOCP F^^ hcocp 
and the above sum. For Fex,Hcocp we use our equation of state (|23|) and equation of state 



for OCP from (this expression for Fex,ocp niay be obtained from Eq.(23) for rj = 0). For 



Fex,hc we use the Carnahan-Starling Eq.(|TB|). In Fig.3 we compare the "complete" equation 
of state (solid lines) with the approximate, based on the assumption of additivity of the 
hard-core and electrostatic parts (dashed lines). As it follows from Fig.3 for small F the 
values of the free energy from the "complete" and approximate equations of state may differ 
significantly; they may even have different signs. The relative difference becomes smaller 
with increasing F (see inset to Fig.3) and decreasing packing fraction. For F > 50 the ap- 
proximation reaches the reasonable accuracy: the error does not exceed 5% for all 77. For 
small packing fractions, 7/ < 0.1 one has the same accuracy already for F > 10. 
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4 Results and Discussion 



A "first-principle" equation of state for tlie classical hard-core one component plasma is 
obtained that has a correct Debye-Huckel behavior at the limit of small plasma parameter F 
and small packing fractions [T —>■ 0, r] 0). It demonstrates a linear dependence on F for 
any packing fraction if F ^ 1. The obtained coefficient 0.881 at the linear leading term in 
this case is close to the corresponding coefficient 0.899 found for the one component plasma 
in the Monte Carlo simulations. The simple analytical expression for the excess internal 
energy reproduces the available MC data with an accuracy of 1 — 3% for the most range of F 
and 7]. The maximal deviation of 12% is observed for small F and rj, where the MC method 
is not very accurate, and we argue that such deviation does not reflect the accuracy of our 
equation. 

To derive the equation of state we apply the Hubbard-Schofield transformation to obtain 
the field theoretical Hamiltonian for the HCOPC and use the Gaussian approximation. The 
Gaussian approximation assumes that one can neglect all terms in the effective Hamiltonian 
which contain the power of field higher than two. As it was shown by the field theoretical 
calculations for the case of the one component plasma (without the hard-core) |T^, the 
corrections to the Gaussian theory are rather small, and we expect that the same is true for 
the system of interest. Physically, this implies that the higher-order coefficients in the field 
theoretical Hamiltonian are smaller than that at the second order of the field. In this case 
the contribution to the field integral from the field-amplitude domain where the higher-order 
terms dominate is small: The main contribution comes from the domain where the Gaussian 
term prevails. Generally, the Gaussian approximation fails near the critical point, where the 
Gaussian term of the effective Hamiltonian becomes small and even vanishes |^^. For the 
HCOCP we do not expect any criticality and thus the Gaussian approximation is expected 
to be valid. 

To obtain the analytical expression for the equation of state we also approximate the 
Fourier transform of the direct correlation function of the reference hard-core system C2(k) 
by its small-/c expansion. Although such approximation deviates from the actual dependence 
of C2(k) at larger k, this does not eventually affect the equation of state: Owing to the long- 
range nature of the Coulombic interactions, only the small-fc behavior of C2(k) is important; 
this is correctly reproduced by the approximation. 

We also analyse the validity of the widely used approximation, where the free energy 
of the HCOCP is represented as a sum of the hard-core and electrostatic component. We 
show that such approximation is rather accurate for small packing fraction and large plasma 
parameter. In the opposite case of large t] and small F the excess free energy may not be 
adequately represented by a sum of a hard-core and electrostatic part. 

Thus we conclude that in general case one has to use the complete equation of state. 
The proposed one possesses a reasonable accuracy for the whole range of parameters for the 
system of interest. 
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Figure 1: Shows the dependence of the excess internal energy of the HCOCP Uex/ NksT 
on the plasma parameter F = Is/dc {Ib = e^/ksT, Oc = (3/47rp)^/^) for different values of 
the packing fraction rj = {'n-/Q)d^p .The curves from top to bottom correspond to rj = 0.001, 
rj = 0.005, rj = 0.020, r] = 0.100 and rj = 0.400, respectively. Points give the Monte-Carlo 



data p9|, |32[: circles correspond to ?7 = 0.001, squares to r] = 0.005, triangles to ?7 = 0.020, 
diamonds to ?7 = 0.100 and down triangles to ?7 = 0.400, respectively. In the inset the same 
dependence is shown for larger range of F. 
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Figure 2: Gives the relative error of the analytical expression (^Sf) for excess internal energy 
of the HCOCP Uex/NkBT as a function of the plasma parameter F. Notations are the same 
as in Fig.l. 
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Figure 3: Compares the "complete" equation of state, Fex,HCOCP (p3|), for the free energy of 
the HCOCP (sohd hues) with the approximate one, given as the sum of the free energy of 
the point particle OCP and that of the hard-core fluid, Fex,ocp + Fhc (dashed lines). The 
curves from top to bottom correspond respectively to r/ = 0.42, rj = 0.2 and rj = 0.05. In 
the inset the same dependence is shown for larger range of F. 
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